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ABSTRACT 

- The use of cellular phones has become widespread. It is predicted that cellular 
phone use will soon become almost as high as normal land-line telephones. The ability to 
track these phones has obvious advantages for intelligence gathering. The problem with 
tracking these phones is that their use is intermittent and they are very low power 
emitters. The use of Unmanned Aerial Vehicles (UAV’s) ould help in the detection of 
these signals. They would extend the capabilities of a ground unit. Once these signals 
have been detected, the best algorithm for tracking them is the Kalman filter because of | 
the intermittent and noisy nature of the received signals. The extended Kalman filter is 
used because of the nonlinearities present in the system. The exploitation of cellular 
signals by using UAV's and the extended Kalman filter is an important framework for 
future use of UAV's in unconventional ways. This thesis explores the use of UAV's to 
exploit seiidias emissions and locate the caller. The time difference of arrival (TDOA) of 
emissions is the fiain method of tracking using the Kalman filter. Further research and 


directions of interest will be proposed. 
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| I. INTRODUCTION 
A. BACKGROUND 

Cellular phones are becoming an everyday appliance to many people. The 
widespread use of these phones has sparked interest in the intelligence community, 
because cellular phone emission characteristics are such that the signals can be tracked. 
By locating a certain phone, the operator's position is likewise known. Being able to 
locate individuals to a certain degree of accuracy could be very beneficial. Currently 
there are few ways to track these phones in a passive way. This work will focus upon the 
problem of finding the position of a cellular phone. 

The UAV (Unmanned Aerial Vehicle) has been very useful in providing 
battlefield surveillance, target information, and naval gunfire spotting. The use of UAV's 
allows the imtelligence field the flexibility to obtain infbemnation that is too risky to obtain 
with manned aircraft or impossible to get with a satellite. The UAV provides some 
measure of stealth since its small size makes it a hard target for radar to detect. The small 
size also is of benefit when the UAV becomes a target, making it harder to hit. 

The UAV could be a valuable asset in tracking cellular phones. A small UAV 
sci in cits alc ol oie chm ell Site, "hs GAN, nition Ho 
| simple omni-directional antenna, would be able to‘ receive cellular transmissions. By 
fitting the sepiaibiia vehicle with a similar auierinal one could measure the time 
difference of arrival (TDOA) for the reception of a cellular emission between the base 


receiver and the UAV. This TDOA is used to estimate the position of the phone. The 





addition of the UAV to the intelligence collection vehicle gives the operator a second 
receiver, which makes the use of these calculations possible. 
B. CELLULAR PHONES 

The excellent work of Lee is used as the main reference when dealing with 
cellular phones [Ref. 1]. Some consider it the bible of cellular phones. 

In general, cellular systems are classified by operating frequencies and channel 
spacing. All systems operate at either 450 MHz or 800 MHz and have channel 
bandwidths of 30, 25, or 20 kHz. There are also two basic types of miielien of voice 
data for these systems, analog and digital. Today, digital systems are just being 
introduced in the commercial market. Digital systems allow a greater number of callers 
a have lower interference problems than tig analog systems. 

_ The system of choice in North America is the Advanced Mobile en Service 
(AMPS), which was first placed into service in 1978. The rest of the world uses various 
systems that are similar to AMPS, but not compatible. This work will focus on gathering 
intelligence from the AMPS cellular system. The wacthnds employed can be modified to 
work with the other cellular systems of the world. | 
C. CELLULAR PHONE EMISSIONS 

Cellular phone emissions are unique in aan ways. The cellular phone system is 
based upon the phone cell and its associated frequency reuse pattern. A cell is a 
geographic area which is assigned certain frequencies of the total spectrum. The size of 
the area depends upon the expected traffic in the cell. An area of operation is broken 


down into a number of these cells with each cell covering a unique portion of the 


2 











operating area, as illustrated in Figure 1. This is an ideal representation, as most cells 


will not have a circular shape and will be of different sizes. 





Figure 1. Cell Structure. 


The cellular structure shown is based upon a seven cell reuse pattern. What this 
means is that each cell of a reuse pattern is assigned a certain segment of the total cellular 
channels available. A seven cell reuse pattern means that the total number of voice 
channels is broken up seven ways and split among the cells. Cells with the same channel 
number use the same voice channels. It is important to separate cells using the same 
frequencies by a distance that will insure that transmissions will not interfere with each 
other. This interference is called cross-talk. To naan cross-talk, the cells are kept 
small and the power outputs are kept low. The AMPS system specifies that the cell 
radius be from 2 km to 20 km. The use of smaller cells allows a greater number of callers 


to be serviced in the same size subscription area. This cell structure causes a problem 
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when trying to track one phone, since there may be a number of phones using the same 
frequencies in a given area, and the phone will switch frequencies as it travels across cell 
boundaries. This switching of frequencies as the phone moves 1s called hand-off. 

The AMPS system is an analog system. This means that voice information 1s 
transmitted using frequency fled Each channel 1s split into a forward voice 
channel (FVC) and a reverse voice channel (RVC). The forward voice channel! transmits. 
from the base station to the mobile station and the reverse voice channel from the mobile 
station to the base station. The forward and reverse channels have a bandwidth of 30 kHz 
and are separated from each other by 45 kHz. In the AMPS system, forward 
transmissions (station to phone) are from 870 to 890 kHz and the reverse transmissions 
(phone to station) are from 825 to 845 kHz as can be seen in Figure 2. The spectrum is 
made up of a total of 666 channels, which are divided into two blocks, A and B, so that 
competing cellular companies can each use half of the total spectrum as delineated by the 
FCC. The channels 313 to 354 are reserved as control channels, and the rest are used for 
voice transmission, giving each carrier a total of 21 setup control channels and 312 voice 
channels. On 26 July 1986, the FCC niin available spectrum by 10 MHz, which 


added 166 more voice channels. The 21 control channels are frequency shift keyed 


(FSK) digital signals. 
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Figure 2. Cellular Channel Spectrum [After Ref. 1]. 


The control shanniels perform a number of tasks, including channel setup. They 
can be access channels, used to setup mobile originated calls, or paging channels, used 
for land originated calls. These signals are digital, and are modulated with frequency 
shift keying (FSK). Digital messages are sent in a Burst ‘mode with a 10 kbps 
transmission rate. There are also two analog control signals, the supervisory audio tone 
(SAT) and the signaling tone (ST). 

The SAT is of specific interest in this thesis. It is a tone that is added to each — 
FVC when an active link is made, and is placed outside of the bandwidth that transmits 
the voice information, but within the channel bandwidth. The mobile phone detects this 
SAT, filters it, and then returns it on the RVC. This was originally designed to allow the 
base station to determine the phase difference when the SAT returned. This phase 


difference could be used to compute the total round trip time from base to mobile and 











back to the base, which could then be used to calculate distance from the — Station to 
the mobile phone. It is not used for distance calculations now, but serves other control 
purposes. Today most systems determine when to hand-off to the next cell based upon 
the signal strength of the RVC or the power difference of the mobile signal received by 
adjacent cell sites. 

D. CELLULAR PHONE EQUIPMENT 

A typical cellular phone area is divided into a large number of cells which use a 
certain frequency reuse scheme. These cells may not be the same size or shape, 
_ depending upon the propagation characteristics of the emissions from the cell site 
antennae and the expected traffic. Each cell site consists of antennae, radios, control 
units, and data terminals. These cell sites are all linked in various ways to the Mobile 
Telephone Switching Office (MTSO). It is here that all the switching is accomplished 
between the seillae network and the land lines; the calls are routed and connected to the 
satonwide telephone network. The MTSO is the heart of the cellular system and 
provides central coordination and sarainisteation for the system. 

Mobile phones have become very small in the last couple of years. The 
transmitting power 1S ielatively low, with the maximum usually 4 watts. The phones do 
not always transmit at the maximum power level, and usually transmit just the minimum 
power required for the base station to receive the signal. This is done to reduce the 
interference Beaween phones. In the future digital phones will become even more 


widespread, as they are better suited for mobile data communications. 








E. PHONE LOCATION PROBLEM DESCRIPTIONS 

Two basic approaches are used to estimate the location of a cellular phone from 
- TDOA information. In the first approach the SAT is used to estimate the total distance 
the signal has traveled from the base station to the phone and then to the receiver. An 
ellipse of estimated location is then calculated. The second approach iss not involve the 
emissions from the base station. In this method the digital Burst emissions from the 
phone are marked in time. The TDOA between two different receivers is then used to 
estimate position. The two approaches may be combined if both forms of information are 
available. 

For both of these methods it will be assumed that one receiver is based out of a 
stationary van (or some other collection vehicle). In addition to the stationary vehicle, the 
user will have at his disposal one UAV that will be able to receive cellular emissions. The 
UAV will be flown a a certain pattern that optimizes the estimation algorithm. 

F. SAT TDOA PROBLEM DESCRIPTION 

The first TDOA filtering problem assumes that a call is in progress between the 
desired phone sca a cell base station. The detection of this condition involves listening in | 
the area sf estimated operation and waiting for the control signals that setup the call and 
assign the channel frequencies. This is easy if the call is originated from land lines, since 

the cellular sien "pages" the cellular phone by broadcasting a setup message from 
| every cell in the subscription area. When the phone responds, the system knows which 
cell to use. Calls originating from a cellular phone are harder to detect. The unique 


identification aaniber for each phone will be the flag that announces that a channel 








assignment is about to be made. Once the channel is known, the SAT from the FVC and 


RVC can be easily tuned in. In this problem the following assumptions are made: 


1. Enough information is obtained from the received signals to filter out 
other emitters. This reduces the multi-emitter problem into a single 
emitter problem. 7 


2. The receivers are referenced to a common time base. This could be 
done by using the clock in the Global Positioning System (GPS) 
system. 


3. GPS will be used in both the land receiver and the UAV. The 
positions of the receivers will be assumed to be without error since the 
accuracy of the GPS system is better than our desired system accuracy. 


4. The receivers each use one omnidirectional antenna. This will keep 
system hardware as simple as possible. 


5. The data link between the UAV and the intelligence van will be 
assumed to be fast enough. 


6. The UAV will fly some set pattern. The efficiencies of these patterns 
will be discussed. A standard pattern will allow automatic flight 
control and relieve the operator from an additional task. 


7. The phone in question will be stationary for the duration of 
measurements. This will eliminate the need to switch frequencies as 
hand-offs occur. 


8. The exact location of the cellular base station in use will be known 
accurately. 


The ranges in question here are limited to the size of one cell. It will be assumed 
that the cell in use has a diameter of 20 km. 
G. BURST TDOA PROBLEM DESCRIPTION 

The Burst TDOA problem estimates the location of the eiditier by measuring the 
TDOA of an emission from the cellular phone between the intelligence van and the UAV. 
The position of the cell base station is not needed and emissions from it are not used. The 


following assumptions are made in developing the problem: 











1. Some sort of DSP method will be used to detect the digital FSK 
control bursts that the phone periodically emits. 


2. The receivers are referenced to a common time base. This could be 
done by using the clock in the Global Positioning System (GPS) 
system. | 


3. GPS will be used in both the land receiver and the UAV. The 
positions of the receivers will be assumed to be without error, since the 
accuracy of the GPS system ts better than the desired system accuracy. 


4. The receivers each use one omnidirectional antenna. This will keep 
system hardware as simple as possible. 


5. The data link between the UAV and the intelligence van will be 
assumed to be fast enough for our purposes. 


6. The phone in question will be stationary for the duration of 
measurements. This will eliminate the need to switch frequencies as 
hand-offs occur. 


7. The UAV will fly some set pattern. The efficiencies of these patterns 
will be discussed. A standard pattern will allow automatic flight 
control and relieve the operator from an additional task. 


8. The exact location of the cellular base station in use will not be known. 

This aii will also be tested for the area of one cell, a 20 km region. It could 
be used for a much larger region when data bursts can be received regardless of which 
cell the phone ts in In theory this is possible, but it must be remembered that there will 
be a great number of phones using the same control channels. A large portion of work, 


which will not be done here, will be required to be able to find the desired ID number. 
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Il. TDOA ESTIMATION 

A. TIME OF ARRIVAL ESTIMATION AND ERROR MODELS 

The Kalman filter is used to estimate the position of the cellular phone from the 
TDOA observations. The Kalman filter is used to estimate the instantaneous state of a 
system perturbed by Gaussian white noise by using measurements that are related to the 
states, but are corrupted by Gaussian white noise. The filter is a statistically optimal 
Sciunitor In this case, the noise of interest will be the measurement noise, faker than the 
process noise. To use the Kalman filter this noise will be assumed to be Gaussian white. 

This assumption can be made because most real noise is in fact close enough to 
this ideal so that the filter algorithms will work. If the actual measurement noise is not 
white, but is instead colored, it can be whitened by using shaping filters. 
B. TIME OF ARRIVAL MEASUREMENT 

| TDOA measurements between the van and the UAV will be used to estimate the 

location of the phone. Since there is no instantaneous physical link between the two 
receivers, it is important that the clocks in each are synchronized precisely. To minimize 
the drift that will occur between any two clocks, they are synchronized to the GPS 
internal clock . The GPS uses this clock in its own position calculations Further checks 
on clock synchronization could be made over the data link by using any of the time check 
algorithms currently in use in computer networks. This assumes that the link between the 
van and UAV is in the form of digital communication. This link could be formed over 
some sort of wireless link, fiber optic line, or even stored for download upon landing the 


UAV. 
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For the SAT TDOA observations the variations encountered in the GPS clock . 
synchronization should be too small to affect the system. The SAT tone is typically at 6 
kHz. This gives a peak to peak time of 166.6 microseconds. The GPS clock 
synchronization error will be a couple of orders of magnitude baiow this. 

The — TDOA observation accuracy depends upon the type of detection 
algorithm used to mark the time of arrival (TOA). The interference may be so great and 
the signal so ee that either receiver or both may miss the signal altogether. This also 
depends upon what part of the cellular emission is exploited. This algorithm could 
trigger the emission of the ID code from the cellular phone, which would relieve the great 
detection burden, but assumes that the phone is engaged in a call. In the future, cellular 
sae will periodically send a registration ID to the cell site so that the MTSO will be 
able to track phone locations. This type of autonomous registration will be assumed to be 
what the Burst TDOA model is receiving. 

C. SAT TIME OF ARRIVAL ESTIMATION 

The SAT is a 6 kHz tone added to the RVC and FVC. Since it is of relatively low 
frequency, it is possible to mark the reception time from the peaks of the signal. As long | 
as the actual time difference is not larger than the peak to peak time, there will be no 
confusion as to which received peaks match the two receivers. The period of the signal is 


166.7 microseconds. Multiplying this by the speed of light yields: 


3x 10° x 166.67 x 10° = 50x 10° (1). 
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Thus, as long as the difference between the two receivers is less than 50 km, this 
signal is used. This does not mean that the range of detection is limited to 50 km, only 
that the separation between the receivers should be less than 50 km. 

D. BURST TIME OF ARRIVAL ESTIMATION 

The Burst TDOA model receives its time of arrival information after the signal 
processing algorithm determines that it has in fact received the desired ID. Since the 
signal is a digital signal that manifests itself as a FSK modulated signal, any one of the 
established computer communication algorithms could be modified for use. Of 
importance is the point where the algorithm triggers the tume measurement. It would be 
preferable to have it trigger on the first received bit. This would require some type of 
‘memory that might be prohibitive. Therefore, the triggering will occur on the last bit. 

As stated earlier, since the signal to noise ratio will be low and the algorithm may 
not ecouaie the signal, the reception of TOA measurements will not occur at a uniform 
rate, but will be modeled as a random process. 

E. SAT AND BURST TDOA 

The TOA for these signals is approximated as a Gaussian random variable. It is 

assumed that each TOA is independent. The TDOA is calculated by using the following | 


equation: 


TDOA=z,.. -Z.., | | (2) 


where Za, = TOA of the emission at the van, and z,,, = TOA of the emission at the 


UAV. 
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The TOA can be expressed as the sum of a true TOA and a Gaussian white noise 
random variable such as 
Z=Z,+V, | (3) 


where z = the observed TOA, z = the true TOA, and v = the noise present in the 
measurement. | 

The mean of the noise is zero and has a variance of V. The expetied value of 
E[VyanVuav] is zero, so the noise components are uncorrelated. The mean of the TDOA 


estimates is 


E[zvan — Za0v] = E[Zovan ~ Zouav }+ E [Yvan J E[eav ]= Zovan — Zona (4) 
It can be seen that the expected value of TDOA is the difference between the true TOA. 


The variance will be 


Cron =E|(TDOA-E [2m ~Zan D EEK. ~Voar) | 7 (5) 
The two noise sequences are independent, 
ELV an Vuav |= 0 ; (6) 
and 
- Bhpo4 = E[Veun + eww [= Youn + Vawv (7) 


Thus, the TDOA is a random variable with a mean value equal to the true TDOA and a 
variance that is equal to the sum of the variances of the TOA measurements. 
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Il. PHONE LOCATION 


The TDOA is the time difference that a signal or series of signals have when they 
are received by receivers that are separated by some distance. By analyzing the geometry 
of the problems the location of the phone is estimated. This location estimate will not be 
a precise point when only one receiver is used, but will be an ellipse or locus of the 
possible positions. The use of two receivers allows the pinpointing of the phone position 
by looking at the intersection of the two possible position loci (ellipses). 

A. USING SAT INFORMATION 

1. Problem Geometry 

In this problem, the position of the cell transmitter is assumed to be known. The 
eeces measures the TDOA between the SAT from the base station on the FVC and the 
SAT — the mobile phone on the RVC. The geometry for a single receiver can be seen 


in Figure 3. 






receiver 


Figure 3. Geometry of SAT Problem. 





The originating SAT travels from the baie station (B) to the phone (P). The 
phone processes the signal and sends it back. The receiver measures the TDOA between 
the emission from the base station and the reply SAT from the phone. In this problem the | 
following values are known: 

d = the distance between the base station and the receiver, _ 

| xb, yb = the position of the base station, 

xr, yr = the position of the receiver, 

Tg = the TOA for the emission from the base station, 

Tp = the TOA for the seatulal from the phone, and 

_c =the speed of transmission, taken to be the speed of light. 
The following values are unknown: 

Xp, yp = the position of the phone, 

Rep = the distance between the base station and the phone, 

Rpg = the distance between the phone and the receiver, 

Ty = the time the base station originates the transmission. 


The time of origination can be calculated as 


Th = 13 -%- | (8) 
The time it takes to travel from the base station to the phone to the receiver, not counting 


processing time at the phone, is 


1-17, =e *En) a (9) 
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Substituting (8) into (9) and solving for the distance 


Rye +R =(To-(To-S)]e=(4-Te +S] | (10) 
which Aaher reduces to 


Rap + Rog = TDOAxct+d. (11) 


Equation (11) describes an ellipse as shown in Figure 3. A formula can be 
constructed to determine the coordinates of the phone based upon the base station 
position, the receiver position, and the TDOA. 

Of more importance in implementing the Kalman filter is the estimation of the 
TDOA based upon the known position of the base station and receiver and the estimated 


position of the phone. Solving (1 1) for TDOA yields: 


I i 2 R 
TDOA = =| \@p- xb) + (9p yb) +(&p- xr} + Gp-yt) - d| / (12) 
where 
d= (xr - xb) + (yr—yby (13) 
2. Loci of Constant TDOA 
The locus of constant TDOA is an ellipse that has its foci at the base station and 


the receiver. The phone can be in an infinite number of possible locations on the ellipse 


that results in the specific TDOA. The problem is to now find the specific location on the 
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locus. One solution is to give the receiver a directional antenna So that a bearing can be 
measured. Using the ellipse locus and the bearing, the position can be estimated. 

As noted before, this problem utilizes two ee that do not have bearing 
capability. By looking at where the loci of each receiver intersect, an estimate of position 
can be made. Notice that having two receivers does not entirely remove the ambiguities, 
since they intersect in two places (see Figure 4). To pinpoint which of these intersections 
is valid, a third receiver could be used. Since the UAV is mobile, each miesarenent 
taken as it moves is considered to be a different measurement from a separate receiver. 


The two receivers then look like countless receivers. 


UAY loci 





van locus 


Figure 4. Loci of Constant TDOA Intersection. 











B. ORTHOGONALITY OF TDOA OBSERVATIONS 

The Kalman filter estimation depends a great deal upon how orthogonal the two 
— loci of positions are. The more orthogonal the intersection, the less error inherent in the 
system. On the other hand, the more parallel the loci are, the worse the estimation will 
be. This inherent error can be shown by calculating error ellipsoids. This is done later in 
this work. 

It is ‘desirable to have the loci orthogonal to each vies Since one of the receivers 
(the UAV) is mobile, this orthogonality is easily obtained through aan 
C. USING ENCODED BURST INFORMATION 

1. Problem Geometry 

This problem does not-consider the location of the base station, since only 
transmissions — the phone will be processed. The TDOA will consist of the time 
difference nee the reception times of a phone’s emissions between two receivers. 
The TDOA is a function of the distance between the receivers. The problem is illustrated 


in Figure 5. _ 
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oe 
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van y 
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x 
Figure 5. Pulse TDOA Geometry. 
The distance from the phone to the van is 
R,.. = V(xp -xvan) +(yp—yvan) . (14) 
Likewise, the distance between the phone and the UAV is 
R,,, = \(xp—xuav) +(yp—yuav) . (15) 
The TDOA formula is therefore 
(16) 


rood = Rem = Raw) 
c 


where: TDOA = the time difference of arrival, 


xp,yp = the phone position, 
20 








Xyans¥ van = the position of the van, 

XuaveYuav = the position of the UAV, and 

c = the speed of light. 

2. Loci of Constant TDOA 

The loci of constant TDOA are hyperbolas with each receiver acting as a focus 
(see Figure 6). There is only one locus that passes through the middle of the two 
receivers. Here again, there are an infinite number of positions where the phone could be. 
In order to resolve this conflict, three receivers would allow the phone to be pinpointed 
by finding where the three loci of position, which pass between each pair, intersect. Only 
two receivers are used in this problem, but the mobility of the UAV allows the simulation 


of more than two receivers by taking readings at different positions for the UAV. 


Loci of constant TDOA 


lic 


Figure 6. Loci for Burst Data. 
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3. Orthogonality of the TDOA Observations 

Kalman filter accuracy depends on how orthogonal the loci of position are at the 
point of intersection. As mentioned salts the more orthogonal, the lower the expected ‘ 
error. This will be shown using error ellipses. 

4. Moving Phone 

In this work, the phone is assumed to be stationary. If the phone moves, a 


dynamic state is assumed. The Kalman filter may still be able to track the phone. 


22 











IV. KALMAN FILTERING 
A. THE EXTENDED KALMAN FILTER 
As will be demonstrated later, the plant equations employed are nonlinear. For 
this reason the normal Kalman filter cannot be used, and the extended Kalman filter must 
be substituted. This involves the use of partial derivatives as linear approximations of the 
nonlinear models. These partial derivatives are evaluated at the estimated values of the 
state variables and are used to calculate the Kalman filter gains. 


For the model the equations of state are 


x(k +1) = £(x(k),w(k),k), | (17) 


where: x(k) = the state of the system, 
w(k) = the plant driving noise, 
f(k) =a function of the states, noise, or k. 


The observation equations are 


~-2(k) = h(x(k),v(k),k), (18) 


where: z(k) = the system measurements, 

v(k) = the measurement noise, 

h(k) = a function of the states, noise, or k. 

Notice that both f(k) and h(k) can be nonlinear pinctioné In both of the problems 
h(k) is nonlinear. The two noises, w(k) and v(k) are assumed to be white noises with zero 


means and covariance's of W and V, respectively. 
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The extended Kalman filter equations are similar to the Kalman equations, as 


shown in the following equations. 


&(k+I|k)=F(%(k\k),k) (19) 
2(k + Ilk) = h(&(k + I]k),k) (20) 
R(k+I k+l) =R(k+1k)+G(k+D[2(k+1)-2(k +k], (21) 


where: X(k|k) = the current state estimation, 

K(k + IIk) = the predicted state estimation, 

Zz(k + ik) = the predicted measurement 

z(k + 1) = the observed ial 

Sik ys the alten vail 

The i is predicted by using (19) and the current estimation of the state. The 
noise is in z(k + 1). Equation (20) predicts the next measurement using the estimate of 
the predicted state and the measurement equation h( ). The final Equation (21) corrects 
the estimate of the state based upon the innovation and the Kalman gains. The innovation 
is the difference between actual cael and the predicted measurements. 

The Kalman gains must be computed on-line in the extended Kalman filter rather 
than off-line as in the ordinary Kalman filter. Here is where the linear approximation of 


the nonlinear model enters into play. First, the Jacobians 
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_ of(x(k), u(k), w(K), k) 
- ax(k) 


o> 





x(k)=<(k;k).u(k )=u(k)w(kj}=0 ? (22) 


dh(x(k), v(k),k) 
ax(k) 


H = 





x(k)=x(k{k),v(k)=0 ° (23) 


The gain equations use these linearized values: 


P(k + Ik) = bP(K{k)O" +Q, (24) 
G(k +1) = P(k + I|k)H™ [AP(k +I|k)HT + R}" (25) 
P(k + Ik +1) = I ~G(k + A P(k +1|k), (26) 


where: P(k|k) = the estimated covariance matrix for the system state following the ," 
measurement update. 


P(k + I|k) = the estimation of the covariance matrix for the predicted state prior 
to measurement update. 


Note that Equations (19) through (26) specify the discrete extended Kalman filter 
- equations. Initial estimates must be provided for the state and the covariance matrix. 


These estimates are usually chosen as 


(00) = E[x(0)}. (27) 


P(0|0) = E [feo — (0) }4(0) - (00) yF Covfx(0)]. (28) 
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It must be noted that since (22) and (23) are linearized about the current estimate 
of the state, the Kalman gains must be computed at the same time observations are being 


made: in real time. 


The Kalman filter guarantees optimal performance, stability, and convergence. 
The extended Kalman filter makes no such promises. It is, in fact, a suboptimal filter. 


The linear approximation that the Jacobian evaluations represent can lead to a divergent 


filter. 
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V. SAT SIMULATION RESULTS WITHOUT DELAY 
A. EXTENDED KALMAN FILTER EQUATIONS 
| The extended Kalman filter will be applied to the SAT problem. The algorithm 


will be tested by simulating it in MATLAB. The following assumptions are made: 


1. The locations of the UAV, the van, and the cell base station are known 
precisely. 


2. A system will be in place to measure the phase differences in the SAT 
from the cell and the phone and convert these into TDOA 
measurements. To this measurement a white Gaussian noise will be 
added. 


3. The phone is not moving. 
4. All locations are referenced to a geostationary coordinate system. 


5. The delay introduced by the phone in returning the SAT will be 
neglected. | 


These assumptions will help to simplify the problem. The delay introduced by the phone 
will be handled later as it is nontrivial. The state and observation equations follow. 


Referring to Figure 3: 


. a, fat k | 
uk +) = Ri) =| ‘| (29) - 


Rop = /(Kt- xc) +(§t-ye) ; (30) 
Ren = y(St—xr) +(ft-yr) . mca) 


d= (xe- xr)? +(yo-yr)* , (32) 
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TDOA = Zp = “IRap +Rog -d]+v, (33) 
Co. 


where: xt, yt = the estimate of the phone's coordinates, 

xr, yr = the receiver coordinates (the van or UAV), 

xc, yc = the cell base station coordinates, 

c = the speed of light, 

v = white measurement noise. 

Notice that since the phone 1s stationary, the state pansion matrix is the identity 
matrix. The nonlinearity of the observation equation (30-33) is what leads to the use of 
- the extended Kalman filter. In implementing the extended Kalman filter, each receiver 
will be treated separately. .The TDOA for the van will be calculated from the current van 
position, and the estimated phone position. This will be used by the filter to update the 
predicted piane position, which is fed into another Kalman filter handles the TDOA for 
the UAV. Each filter works off the results of the prior filter. In this way, the two stage 
filter should converge estimates upon the true position of the phone. A two stage filter 
| algorithm is used, rather than trying to combine both into a single stage filter 
To use the extended Kalman algorithm, the Jacobian of h( ) must be found. 


Taking the derivatives and making substitutions, 








nx=_ Xt-Xe | XtoXT XI XC 7 (34) 
cl Rup Ror d 
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l| yt- yt— ~ 
pyod Soe Je we 


Ruy “Re «© jie ae 
H=[hxhy]. (36) 
Equation (36) can then be used in the extended Kalman filter equations 
P(k+ilk & P(k|k}#Q, (37) 
P(k+I[k+1}- [I-G(k+1)H(k) P(k+|k), (38) 


G(k+1)=P(k+I[k Ak) [AQ )P(K+IKAQ)S +R YP (39) 


Note: A(k) whould be H(k+ I|k if the phone was not stationary. 


For stability, the Joseph form of the covariance update equation will be used 


instead of Equation (38). 


B= [i -~G(k + NA(K) |, (40) 


P(k+I|k+1)=BP(k+1|k)B* +G(k+1)RG(k+1)". (41) 


The use of the Joseph form is more computationally expensive, but is less 


sensitive to round-off errors and will not lead to negative eigenvalues. 
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B. TDOA SIMULATIONS 

To test the algorithm the following scenario was used. The cell base station is 
placed at the coordinates (10,10). The intelligence van is placed at (10,0). The UAV is 
free to fly a circular pattern that will be centered either on the cell center or the van. The 
phone is free to be placed inside or were the pattern of the UAV. The initial error’ 
covariance matrix and the variance of plant excitation are varied until the filter reaches its 
best convergence characteristics. Arbitrary white noise has been added to the TDOA 
measurements. The initial estimation of the phone position will always be (8,8). There is 
no particular reason for this and the initial guess is not plotted. 

For each scenario three plots are shown. The first shows a time vs. state 
estimation for the phone. This shows the convergence characteristics of the filter. The 
second plot will be an X-Y plot showing the trajectory of the estimates of phone position. 
This plot will also show the estimated loci of position for every fifth observation. This 
gives an idea of i orthogonal the loci are. The UAV will be represented by an 
asterisk, the track of the estimated position by a jagged solid line, and the loci of position 
by dashed lines or ellipses. The ellipses of three sigma will be shown as a solid ellipse. 
The third plot will be a close up of the area around the actual phone position to 
investigate how accurate the filter is. 

The observations were first taken at a uniform interval on the UAV trajectory. 
The results of this were less than satisfying. The measurements were then taken at 
random intervals of the UAV trajectory, which worked much better. The MATLAB 
programs used are included in the appendices. 
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1. Scenario One 
In this scenario, the UAV is allowed to fly a circle around the cell center with a 
radius of 10 km. The phone is first placed at (13,14). The best convergence occurred 


with 


5 0 2 5 0 i0 20 
oe 7 x107, PO= - «10°, andR=1.694x10°. (42) 


It is evident in Figure 7 that the filter saavetaes quickly--within 8 observations. 
After convergence, the filter remains stable. 

In Figure 8 it can be seen that the trajectory of the estimation of the phone 
position does a fairly good job at converging on the phone location. It does show some 
large perturbances, but settles down on the phone location quickly. Figure 9 shows that 
the filter does wander around the location of the phone a little. The final two three sigma 
error ellipses are shown in the close-up. These show that there is a probability of 98% 
that the phone is in this ellipse. The final ellipse has axis of dimension 180 m by 40 m. 


This is very acceptable. 
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Estimate of the x and y coordinates 
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_ Figure 7. Observations vs. Phone Location, Scenario One. 
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Figure 8. Trajectory of the Phone Estimation, Scenario One. 
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Figure 9. Close-up of Trajectory, Scenario One. 
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2. Scenario Two 

The next test keeps everything the same except the phone position. The phone is 
now placed at (10,24), outside of the UAV pattern but in-line with the cell center and the 
van. The results are shown in Figures 10, 11, and 12. 

The final accuracy of the filter is not as good as in Scenario One. It takes the 
filter just over 20 observations to converge. Figure 11 shows how the ellipses of location 
intercept each other around the phone location. The intersection is not as orthogonal as 
before, which would cause the longer ecaversence time. The final error ellipse is slightly 


larger than before, at 180 m by 50 m. The filter is rather well behaved. 
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Estimate of the x and y coordinates 
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Figure 10. Observations vs. Location, Scenario Two. 
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Figure 12. Close-up of Trajectory, Scenario Two. 
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3. Scenario Three 

As a final test the phone is placed at (25,11). The filter converges within 5 
observations, as seen in Figure 13. Figure 14 shows that the intersection of the ellipses of 
location are fairly orthogonal. This is more visible in the close-up of Figure 15. The 
major axis of the ellipse is aligned parallel to the location leas which is mpetted 
The size of the error ellipse is 30 m by 350 m. Again, this is acceptable. It should be 
noted that even though the ellipse has such a large major axis, the actual track stays well 
within the center of the ellipse and over the phone's position. If only this small section of 
the on is inspected, the accuracy is quite good. Obviously the geometry of the problem 


results in loci of position that are not as orthogonal as the two other situations. 
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Figure 13. Observations vs. Location, Scenario Three. 
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4. Results 
The algorithm worked fairly well. The best accuracy occurred when the phone 
was placed within the UAV pattern. The worst convergence occurred when the phone 
was placed outside the UAV pattern and not in-line with the van and cell base. The 
ellipses of error show the relative nature of the final accuracy’s for the filter. 

Observations were made at random intervals on the UAV track, which helped the — 
filter to converge both faster and more accurately. Although not shown, when 
observations were taken at fixed intervals of its flight path, the filter had much worse 
performance. One interesting result of using random observations is that each simulation 
run on the same initial conditions gave different results. The plots shown are 


representative of the majority of these runs. 
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VI. SAT SIMULATION RESULTS WITH DELAY 

A. EXTENDED KALMAN FILTER EQUATIONS 

In the last chapter the extended Kalman filter was applied to the SAT problem 
without consideration to the delay encountered at the phone. Each phone has some delay 
associated with processing signals. This is usually set by the manufacturer. If the delay 
was known, it would be a simple matter to add it into our TDOA equations. It must be 
assumed that it is not known, and as an unknown variable it can introduce unacceptable 
error into our filter. Figure 16 shows the trajectory that results when using the filter 
developed in Chapter V. A delay of 4 microseconds was arbitrarily chosen. 

The position of the phone was (13,14), which previously yielded good results. 
The filter does not convetge upon the actual phone position, but is offset by the error 
introduced by the delay. Notice that the y coordinate of the filter converges to almost 15 
km. The daha view in Figure 17 further shows how the delay affects the filter result 


It does not converge around the phone as before. 
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Figure 16. Trajectory with a 4 Microsecond Delay. 
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Figure 17. Close Up of Uncompensated Delay. 
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To compensate for the unknown delay, the filter equations must be modified. The 
unknown delay is now considered a third state of the equation, and is to be estimated. 
This is accomplished rather easily. The new state is | 

xt(k) 
x(k+1)=x(k)=Hyt(k)}. | (43) 
d(k) 
The TDOA equations are modified by adding d to the real observations and d(k) 


to the estimated TDOA. The Jacobian 1s also modified by 


H=[hxhy ht], oo (44) 


where hx and hy are the same as before. | 


1. Scenario One 
The initial estimate of delay will be 5 microseconds. The initial R stays the same 
as before. The challenge with this algorithm is to find the Q and PO that give the best 


- convergence. After much experimentation the following was used. 


5x10° 0 0 5x10" 0 0 


Q=| 0 5x10 0 and PO=; O 5x10" 0 (45) 
0 0 2x10 | 0 0 25x10"! 


The results of this filter are seen in Figures 18, 19, and 20 for a phone position of 
(13,14) and an actual delay of 4 microseconds. Figure 18 shows the convergence of the 
states within 10 observations. In Figure 19 it can be seen that the error ellipsoids start out 


very large and become smaller as the filter converges. The close-up in Figure 20 shows 
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the final error ellipsoid has the dimensions 1,300 m by 80 m. This is larger than the 
results shown in Chapter V. The filter converges upon the phone's position rather well. 


[t is much better than the results with no compensation for the delay. 
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Figure 18. State Plots with 4 Microsecond Delay, Scenario One. 
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Figure 19. Trajectory for 4 Microsecond Delay, Scenario One. 
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Figure 20. Close-up of Trajectory for 4 Microsecond Delay, Scenario One. 
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2. Scenario Two 

Here the phone 1s placed at the position (10,24) as was done in the previous 
chapter. The results are seen in Figures 21 and 22. 

In this sample. the delay converges to 25 microseconds and throws off the 
position estimates. The y coordinate of the estimated phone position converges to 20 km, 
4 km off the mark. The state trajectories in Figure 22 show just what effect the error in 


delay estimation has on the position estimation. 
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Figure 21. State Plots with 4 Mircosecond Delay, Scenario Two. 
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3. Scenario Three 
The phone will be placed at (25,11) which gave the worst results for the filter 
discussed in Chapter V. The results are seen in Figures 23, 24, and 25. The filter 
converses relatively well, coming to a steady state within 10 sieenations The delay 
state is still not converging to the actual state, but 1s much closer than before. 
The major axis of the — ellipsoid is large at over 2 km. Note, however, that the 
actual track converges to a point around 20 m to the west of the phone, as seen in Figure 


25. If one was looking at the track alone, this filter would give acceptable results. 
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Figure 23. State Plots with 4 Microsecond Delay, Scenario Three. 
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Figure 24. Trajectory for 4 Microsecond delay, Scenario Three. 
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Figure 25. Close-up of Trajectory for 4 Microsecond Delay, Scenario Three. 
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4. Results 


The above plots were chosen because they were the best results for a number of 
runs Mitvihe same initial eondidons. Each run gave slightly different results. This filter 
gave a number of results that did not converge, and thus the filter depends upon the 
random measurements for its stability. While the number of results that did not converte | 
was not great the filter employed in Chapter V always converged. 

| By treating the dcisy as a state to be estimated, we have solved the problem of 
most of the offset introduced by an unknown delay. The filter did not always converge 
on the actual delay. It usually got close enough so that the introduced offset error was 
acceptable. More adjustment of the Q and PO matrices may make the delay estimation 


always converge to the real value. 
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Vil. FSK BURST TDOA SIMULATION RESULTS 
A. EXTENDED KALMAN FILTER EQUATIONS 
The extended Kalman filter is applied to the Burst TDOA problem. The 
algorithm will be tested in MATLAB for convergence and stability. The following 


assumptions are made: 


1. The locations of the UAV and the van are known precisely (GPS in 
each). The location of the cell antenna is not needed. 


2. A signal processing system will be in place to filter out the desired 
phone registration number and measure the arrival time to compute the 
TDOA measurements. To this measurement a white Gaussian noise 
measurement will be added. 


3. The phone is not moving. 


4. All locations are referenced to a geostationary coordinate system. — 


These assumptions help to simplify the problem. The estimated state equation is 


the same as before. The predicted observation equation is now 


-2(k + Ik) = ~ (Rov 2h), (46) 


where: R,,,, = distance between estimated phone position and the van, and 

Reavy = distance between estimated phone position and the UAV. 

Here the observed TDOA is the time that it takes the coded burst from the phone 
to reach the van and UAV. The Kalman equations are the same ms those employed 


before. The Jacobian has changed slightly, to 
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1|}{ Xt-xvan xt-xuav){ yt-yvan yt- | 
H(k) = — tava = Stet) Ho = Pe ee oad pe (47) 
. Ravan R pvav 


The Joseph form of the covariance update is again used to minimize round off 
errors. The filter algorithm works the same as the SAT algorithm did. The difference 
here is that only one loci of estimated position can be calculated for each TDOA 
measurement. The filter only has one stage, which processes the single observation. It 
then loops around for the next measurement. It is noted that if the UAV did not move, 
the loci of estimated position would not move, and no clear interception would be found. 
The loci would also not be very orthogonal, and a great deal of error would be present if 
it did eonvense . 

Like before the Q and PO matrices are varied until the filter gives the best 
convergence. The van and UAV locations are as before. The phone will be placed in 
‘three pattie as eles 
B. TDOA SIMULATIONS 

In this case, all of the initial conditions will be the same as before. The best 


convergence occurred with 


10° 0 5x10" 0 a 
= and PQ= . . 48 
Q F 0 sae Sa 


The simulation is allowed to run for 40 observations. The observations occur at 
random intervals of the UAV track. This not only gives better convergence, but simulates 


the noisy environment and difficult detection of these bursts. The three scenarios are run 
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as before with the same resulting three plots shown. The lines represent the loci of 
estimated position in the vicinity of the estimation. These will show how orthogonal the 
measurements are. 

1. Scenario One 

The phone is placed inside the UAV pattern at (13,14). The results are seen in 

Figures 26, 27 and 28. | 

This filter acted very stable. It is seen from Figure 26 that the filter converges | 
within 30 observations. The trajectory shown in Figure 27 een a no nonsense track 
and heads directly for the phone's position. This is seen clearer in Figure 28 as the track 
converges directly on the phone's location. Although this algorithm is slower to 


converge, it appears to be more accurate than the SAT algorithm. 
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Figure 26. Burst TDOA Problem, Scenario One. 
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Figure 27. Trajectory of Burst TDOA Problem, Scenario One. 
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Figure 28. Close-up of Trajectory, Scenario One. 
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2. Scenario Two 

In this scenario.the phone is placed at (10,24), outside of the UAV sateu but 
opposite the van. The results are as shown in Figures 29, 30, and 31. Figure 29 shows 
that the convergence of this situation is much slower than before. It can be seen that the 
filter in fact has not fully converged by the 60" observation. Figure 30 shows the track 
heading in the general direction of the phone, getting ever closer. The close-up in Figure 
31 shows that the track approaches to within 10 m of the phone's position. This is very 
accurate. 

The filter seems to try to converge but never makes it. This phone position also 


gave a great number of results which did not converge at all. 
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Figure 29. Burst TDOA Problem, 


Scenario Two. 
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Figure 30. Trajectory of Burst TDOA Problem, Scenario Two. 
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Figure 31. Close-up of Trajectory, Scenario Two. 
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3. Scenario Three 

The final test places the phone outside of the UAV pattern at (25,11). The results 
are shown in Figures 32, 33, and 34. In general, the new position gives the same results 
as found in Scenario Two. The change does not have much of an effect upon the filter. 

This filter never converges around the exact phone location. It displays the same 


single mindedness by constantly getting closer to the phone. 
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Figure 32. Burst TDOA Problem 


, Scenario Three. 
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Figure 33. Trajectory of Burst TDOA Problem, Scenario Three. 
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Figure 34. Close-up of Trajectory, Scenario Three. 
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4. Results 


The algorithm worked quite well when the phone was placed inside the UAV 


track (Scenario One). The random nature of the observations again caused differing 
results with the aie initial conditions. This algorithm has the advantage that the cell site 
does not have to be known, and only one frequency need be monitored. The results are 
close to the actual phone position. In fact, when the phone was placed inside the UAV's 
‘flight circle, the results were very accurate. This method favors having the UAV fly 


around the suspected position. 
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VII. COMBINED SAT AND BURST FILTER 

A. ALGORITHM CHANGES 

In this chapter the two mieitiods of finding the position of the phone will be 
combined into one filter. ils alate that both the SAT information and the Burst 
information are available, and that the location of the cell station is known. This method 
is easily implemented by taking the last two filters and combining them n series. 
Basically, in one loop of the program, the filter will first modify its position estimate 
based upon ihe SAT information. It will then feed this estimate to the filter produced for 
the Burst filter. The algorithm then loops back to the beginning using the Burst updated 
estimation. This is certainly the easiest implementation of this idea. The program could 
also have been modified to use both measurements of TOA at the same ane To keep 
things simple, it is assumed that the SAT delay is ‘is for the SAT calculations. This has 
been ene to test if these two algorithms can work together on a basic level. 
B. SIMULATIONS 

The simulations will follow the same pattern used in the previous chapters. The 
same locations are used so that the results may be compared to prior results. Randoni: 
mieasurenicals will again be assumed on the track of the UAV. This is not so unrealistic, 
since in a real scenario the receivers may be able to discern only occasional pieces of 
information from the cellular phone. As the phone travels through changing landscape, 
each receiver may temporarily find itself in a blind spot of the signal. This algorithm will 


test how well one can predict the phone position when the SAT and Burst information are 
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spotty. Hopefully, the combined information will be enough to let the combined filter 
converge. 

1. Scenario One 

The phone is placed :at (13,14) to test how the algorithm reacts to having the 
phone inside of the UAV pattern. The results can be seen in Figures 35, 36, and 37. In 
Figure 35 it is seen that the filter converges within 4 observations. This is very quick in 
relation to the other filters. Figure 36 shows how the filter heads straight for the phone's 
location. The close-up in Figure 37 gives a better idea of what is happening around the 
point of convergence. The three sigma ellipses are 160 m by 100 m. This is not 
significantly smaller than the other results, but when the plot is examined it is evident that 
the filter stays within 20 m of the phone. This sort of accuracy is what is needed for — 


practical application. 
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Figure 35. Burst TDOA Problem with Combined Filter, Scenario One. 
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Figure 36. Trajectory of Burst TDOA Problem with Combined Filter, Scenario One. 


80 





(km) 








39 12.92 1294 1296 12.98 13 13.02 13.04 13.06 13.08 13.1 
(km) 


Figure 37. Close-up of Trajectory with Combined F ilter, Scenario One. 
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2. Scenario Two 

As in the previous chapters, the phone is placed at (10,24) to test how the 
algorithm handles the situation with the phone outside of the UAV track and in-line with 
the cell center and the van. 

Figure 38 shows the filter converging within 15 observations. It is seen in Figure 
39 that the filter makes the largest jumps within the first 4 observations. Figure 40, the 
close-up of the phone, shows the three sigma ellipse to be 200 m by 70 m. The final 
convergence does not leave the phone's position by more than 40 m once it has settled 


down. 
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Figure 38. Burst TDOA Problem with Combined Filter, Scenario Two. 
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Figure 39. Trajectory of Burst TDOA Problem with Combined Filter, Scenario Two 
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Figure 40. Close-up of Trajectory with Combined Filter, Scenario Two. 
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3. Scenario Three 

In this scenario the phone is placed at (25,11) to test the response of having the 
phone outside the UAV's track and not in line with the van and cell station. 

Figure Al shows that this filter converges quickly in this scenario, at less than 10 
observations. The track of the filter shown in Figure 42 shows that the filter has some 
trouble in the first couple of observations, estimating the position to be as far as 65 km 
away from the actual position. What is also interesting is that the estimates stay outside 
of the UAV track. The close-up shown in Figure 43 shows the three sigma ellipse having 
the dimensions of 160 m by 120 m. Once converged the algorithm keeps its estimates 


within 40 m of the phone's actual position. 
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Figure 41. Burst TDOA Problem with Combined Filter, Scenario Three. 
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Figure 42. Trajectory of Burst TDOA Problem with Combined Filter, Scenario Three 
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Figure 43. Close-up of Trajectory with Combined Filter, Scenario Three. 
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C. RESULTS 

The biggest benefit of combining the SAT and Burst filters is not shown in the 
plotted figures, as the combined filter eae every time it was run. The stability of 
this filter far surpasses the stability of its two components working separately. Out of the 
50 runs made with this algorithm, only bie failed to converge. Remember that the sein 
of the reading is taken at a random spot on the UAV track so that each run is running off 
of essentially unique data each time. The robustness of this algorithm is of great value. 
The accuracy and time to converge is also as good if not better than the other two filters. 
It converges as fast as the SAT filter and has as good accuracy as either filter. 

It should be noted that both an SAT and Burst data reading was taken at each 
random measurement point. In the future, each point may be limited as to which signal 
can be received. The filter would have to be modified to operate only when dain was 
received, ignoritig the missing data in the loop. The delay estimation would also have to 


be included in the next iteration of this program, as the nonstandard delay of the SAT 


signal can affect accuracy. 
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IX. CONCLUSIONS AND RECOMMENDATIONS 

A. CONCLUSIONS 

All of the extended Kalman filter algorithms developed here perform well. The 
position of a cellular phone can be estimated to a great degree of accuracy by measuring 
the TDOA for certain cellular emissions. The use of a base station and a UAV has shown 
to be highly beneficial to these algorithms. The movement of the UAV causes the 
shifting of the loci of positions for the filters. This shifting of the loci allows the Kalman 
filter to estimate the intersection and thus the phone position. | 

The filters developed are highly dependent upon the type of emission being 
exploited. Each type of emission has strengths and weaknesses with regard to locating | 
phones. The SAT filter is very accurate, and has a good convergence time. The SAT 
signal is continuous and easy to intercept once the channels being used are known. The 
major detection is that the phone has an inherent delay which adds an error constant. 
This thesis has shown that it is possible to predict this delay by using a three state 
Kalman filter. This filter worked fairly well, but was the most unstable filter of the 
group. Many times it would diverge on its estimates or converge to the wrong delay 
value, which does not solve the delay problem. The filter does show that it can work 
well, though, and has the potential to become a very stable and accurate filter, given 
further development. 

The Burst Kalman filter relies on the reception of control codes from the phone. 
This reception and isolation of the signal is a thesis problem in itself. Assuming that this 


information can be obtained, the filter showed that the concept of measuring only the 
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TDOA does work well. This filter showed good stability, but was slow to converge and 
often never actually reached the phone position. 

When the two filters were combined, ignoring the delay of the SAT problem, the 
resulting filter had the best features of both of its components. The SAT Burst filter was 
accurate, quick to converge, and very stable. The results of each filter enhances the 
performance of the next. This combination shows the most promise since one or the 
other signal or both may be available at any one point in time. 

These algorithms have the advantage that the receivers are very simple and 
omnidirectional. The UAV is perfect for such use, since the hardware required by such a 
systee would easily fit. (Though a small plane or boat could also be used). The 
algorithms are not too complex and would not overpower a small mobile computer. 

B. RECOMMENDATIONS _ 

This thests is just the beginning of an application which has great potential. 
Although sul one UAV pattern was investigated here, different patterns may have 
differing efficiencies. A development of these patterns and their uses would have merit. 
The study of how well these filters track moving targets would be a natural progression. 
The a priori estimates have a great effect upon the filter's behavior. A study of how these | 
estimates affect stability, convergence, and accuracy would be sonnweie 

It would be beneficial to develop an algorithm that would be able to set the a 
priori estimates based upon expected conditions and the first observation. A neural net 


might be trained to do this so that the filter gives optimum results. 
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Two receivers could have been used on each vehicle. This might increase 
accuracy. Two UAV's might also be used, and their interaction would be an interesting 
area of investigation. 

The combined SAT Burst filter should be modified to account for SAT delay at 
the cellular phone aad to operate on available data. Some data may be suppressed or 
erroneous. The algorithm should be able to distinguish which is useable and apply it to 
the Kalman filter. 

When all is said and done this idea and technique show merits for the tracking of 
cellular phones. The widespread use of these phones makes this tracking even more 


_ attractive. 
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APPENDIX A. ERROR ELLIPSES 


In Ref. 3, a theory is developed that derives the error ellipsoids from a given 
covariance matrix of estimation error. This is P in this work. The source states that the | 
eigenvalues of P are @,,Q,,...,@, and the corresponding eigenvectors are 


n 


(1) 


ie a 


jy’. The quadratic form of the hyperellipsoid is developed as 


X24 pte np (1A) 


This forms an ellipsoid in n dimensional space, where n corresponds to the 
number of states our filter has and y’ is the linear coefficient of y) and so on. If n is 
greater than two hyperellipsoids are formed. These hyperellipsoids have surfaces which 
have equal probability densities. Thus, for a certain 2 the probability that a point will lie 
inside the ellipsoid can be computed. 


The probabilities have been calculated for a few values of dimension, n and 2: 


Rd 


l 2 5 


we 0.683 0.955 0.997 


0.394 0.862 0.989 
7 3 0.739 0.971 
Table Al. Probability of n Dimension Ellipsoids. 
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Since it is desired to show the accuracy potential of the filters and the dimension 
is n=2 for plotting purposes, 1=3 or a three sigma error ellipse is used. This final 
ellipse yields an ellipse for which there 1s a 99.7% probability that the desired value is 
inside. 

The MATLAB program which does this uses the MATLAB "eig" function to — 
calculate the eigenvectors and eigenvalues of the covariance matrix. It — forms the 
ellipse from Equation 1A and converts this into a polar, which is easier to get locations 


for plotting. 
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APPENDIX B. MATLAB PROGRAMS 


TDOA -- THE BURST FILTER 
% Michael P. Fallon 


% TDOA: This program simulates the TDOA burst estimation problem. 
% The van and UAV are assumed to have single antennas in 
% this scenario. The UAV is allowed to fly a racetrack pattern 
% The reception of the identification signals from the phone are 
% assumed to occur at random intervals on the UAV track. 

% 

% 


clear 


c=3.00e+08; % speed of light 
R1=1.694e-20; ~% The covariance of measurement error 


Q=10e+04*eye(2); % Variance of Plant excitation noise 
PO=5e+20*eye(2); % Initial error covariance 


% Choose to calculate the error covariance of 3 sigma 
C=3.0; | 


tf=60; 
t=(0:tf); 


% the time vector 


len=length(t); 
l=ones(1,len); 


% Initialize storage for Kalman gains and the states 
g=zeros(2,len); 
X=zeros(2,len); 


% Emitter location input by the user. 
xt=input(’ Enter the x coordinate of the phone in km..’); 
yt=input(’ Enter the y coordinate of the phone in km..’); 


xt=xt* 1000; 
yt=yt* 1000; 


% A priori estimate of the location of the target 
X0=[10000; 10000]; 


% Location of the van and UAV 
97 


xvan=10000"!+0"t: 
yvan=0"1+0"t; 


% The UAV will move in a circular pattern centered on the cell. 
dtheta=pi*rand(1,len)/5; | 
theta(1)=-pi/2; 


for i=2:len 
theta(i)=theta(i-1)+dtheta(i); 
end 


% Set the center of rotation and the radius for the UAV. 
Cx=10000; 

Cy=1 0000; 

r=8000; 


xuav=Cx"l+(r*l).*cos(theta); 
yuav=Cy’l+(r*l).*sin(theta); 


% Calculate the TDOA observations for the two receivers 
Z=(1/c)*(sqrt((xt-xvan).42+(yt-yvan).42)-... 
| sqrt((xt-xuav).%2+(yt-yuav).%2)); 


% Inject zero mean noise into the observations 
n=sqrt(R1)*randn(1,len); 
2n=Z+Nn; 


% Initialize the Kalman Filter 
X(:,1)=X0;  % A priori estimate of the emitter location 
PK=P0; % Initial error covariance matrix 


% calculate the filtered estimates of the emitter location. 

for k=1:len-1, | 
Rvan=sart((X(1,k)-xvan(k)).42+(X(2,k)-yvan(k)).42); 
Ruav=sart((X(1,k)-xuav(k)).42+(X(2,k)-yuav(k)).42); 


hx=(1/c)*((X(1,k)-xvan(k))./Rvan-(X(1,k)-xuav(k))./Ruav); 
hy=(1/c)*((X(2,k)-yvan(k))./Rvan-(X(2,k)-yuav(k))./Ruav); 


HK=[hx hy]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 
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Zhat=(1/c)*(sart((X(1,k)-xvan(k))42+(X(2,k)-yvan(k))42)-... 
sqrt((X(1,k)-xuav(k))42+(X(2,k)-yuav(k))2)): 
Zh(k)=Zhat; 


% The new error covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK’*inv(HK*PK*HK’+R1); 


% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; ~ 


% Joseph form of covariance update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)'+GK*R1*GK’; 


% Calculate the smoothed estimate of the bearing to the phone 
X(:,k)=X(:,k)+GK*(Zn(k)-Zhat); 
X(:,k+1)=X(:,k); 


% save the error covariance matrices 
if k==1 | 
P=PK; 

else 

P=[P;PK]; 

end; | 


% save the Kalman gains 
g(:,k)=GK; 


end; 


% Plot the output 


% Plot the estimates of the phone location. 
Figure(1) | 

clg; 

plot(t,X(1,:)/1000,t,X(2,:)/1000) 
title(’Estimate of the x and y coordinates’); 
xlabel(’/(Number of observations’) 
ylabel(’Phone location in km’) 

grid 

pause 
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% Steady state estimation of emitter location 
xtavg=mean(X(1,k-10:k)); 
ytavg=mean(X(2,k-10:k)); 


% Plot the grid 

Figure(2) 

plot(xt/1000,yt/1000,’"*’) 

hold on | 
plot(xvan(k)/1000,yvan(k)/1000,’r*’) 
plot(X0(1)/1000,X0(2)/1000,’+’) 


% Plot the error ellipses. 

for k=1:len-1; 
plot(xuav(k)/1000,yuav(k)/1000,’g*’); 
P1=[P(2*k-1,1) P(2*k-1,2);P(2*k,1) P(2*k,2)] 
[xg, yg,Emax]=ellipa(P1,C,X(1,k),X(2,k)); 
%plot(xg/1000,yg/1000,’b-’); 

end 


plot(X(1,:)./1000,X(2,:)./1000,'r); 
axis(’equal’) 

hold off 

title(’Estimate of phone’) 
xlabel(’(km)’);ylabel(’(km)’); 


pause 


% Plot close up of the estimate of the emitter location 
Figure(3) 

clg | 

axis([(xt-2000)/1000 (xt+2000)/1000 (yt-2000)/1000 ... 
(yt+2000)/1000)); 

hold on 

plot(xt/1000,yt/1000,’*’) 
plot(X(1,:)./1000,X(2,:)./1000); 


k=len-2; 

p=P(2*k-1:2*k,:); 

[xg, yg,Emax]=ellipa(p,C,xtavg, ytavg); 
%plot(xg/1000,yg/1000,’b-’); 


Ixg,yg,Emax]=ellipa(P1 ,C,xtavg, ytavg); 
%plot(xg/1000,yg/1000,’b-’); 
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hold off 
title’Close up of position of the Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 


TDOA2-SAT FILTER WITHOUT DELAY 


% Michael P. Fallon 


% TDOA2 This program simulates the TDOA estimation problem. 
% In this scenario we assume that the van and UAV are using 
% the SAT to calculate the ellipses of location. 


clear 
c=3.00e+08; 
% speed of light 


R1=1.694e-15; % The covariance of measurement error 
Q=5e+02*eye(2); % Variance of Plant excitation noise 
PO=5e+10*eye(2); % Initial error covariance 


% Choose to calculate the error covariance of 3 sigma 
C=3; 


tf=50; 
t=(0:tf); 


% the time vector 
len=length(t); 
l=ones(1,len); 


% Initialize storage for Kalman gains and the states 
g=zeros(2,len); 
=zeros(2,len); 


% Emitter location : 

=input(’ Enter the x coordinate of the phone in km..’); 
yt=input(’ Enter the y coordinate of the phone in km..’); 
d=input(’ Enter the phone. delay in microseconds..’); 
d=d*1e-6; 


xt=xt*1000; 
yt=yt*1000; 
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% A priori estimate of the location of the target 
X0=[8000;8000); 


% Location of the cell site antennae 


— xce=10000'1; 
yc=10000"1; 


% Location of the van and UAV 
xvan=100007!+0*t; 
yvan=0"1+0"t; 


dtheta=pi*rand(1,len)/2; 
%dtheta=pi*.1; 


theta(1)=-pi/2; 


for i=2:len | 
theta(i)=theta(i-1)+dtheta(i); 
end 


Cx=10000; 
Cy=10000; 
r=10000; 


xuav=Cx*l+(r'l).*cos(theta): 

yuav=Cyl+(r*l)."sin(theta); 

%xuav=15*l+1*t: 

%yuav=0'l; 

%o Calculate the TDOA observations for the two receivers 


Z1=(1/c)*(sart((xt-xc).42+(yt-yc).42)... 
+sqrt((xt-xvan).A2+(yt-yvan).42)... 
-sqrt((xvan-xc).42+(yvan-yc).%2))+d; 

22=(1/c)*(sqrt((xt-xc).\2+(yt-yc).A2)... 
+sqrt((xt-xuav).A2+(yt-yuav).42)... 
-sqrt((xuav-xc).42+(yuav-yc).%2))+d; | 

% Inject zero mean noise into the observations 

ni=sart(R1)*randn(1,len); 
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n2=sqrt(R1)*randn(1,len); 


Z1in=Z1+n1; 
Z2n=22+n2: 


% Initialize the Kalman Filter 
X(:,1)=X0; %A priori estimate of the emitter location 
PK=P0; . % Initial error covariance matrix 


% calculate the filtered estimates of the emitter location 
for k=1:len-1, 
% perform the first calculation for the measurement at the van. 


Rtc=sart((X(1 sk)-xc(k)).42+(X(2,k)-ye(k)).A2); 
Rtr=sqrt((X(1 ,k)-xvan(k)).42+(X(2,k)-yvan(k)).A2); 
Rre=sqrt((xvan(k)-xc(k)).42+(yvan(k)-ye(k)).42); 


hix=(1/c)*((X(1,k)-xc(k))./Rtc+(X(1,k)-xvan(k))./Rtr-... 
(xvan(k)-xc(k))./Rrc); 

hi y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yvan(k))./Rtr-... 
(yvan(k)-yc(k))./Rrc); : 


HK=[h1x h1y]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 


Zhat1 =(1/c)*(sqrt((X(1 ,k)-xc(k)).42+(X(2,k)-ye(k)).42)... 
+sqrt((X(1,k)-xvan(k)).42+(X(2,k)-yvan(k)).42)... 
-sqrt((xvan(k)-xc(k)).42+(yvan(k)-yc(k)).42)); 


% The new error covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK"*inv(HK*PK*HK’+R1); 


% Calculate the updated error covariance matrix | 
%PK=(eye(2)-GK*HK)*PK; 


% Joseph form of covariance update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*R1*GK’; 
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% Calculate the smoothed estimate of the bearing to the phone 
X(:,K)=X(:,k)+GK*(Z1n(k)-Zhat1 ); 
X(:,kK)=X(:,k); | 


%e perform the calculations for the measurement at the UAV 


Rtc=sart((X(1,k)-xc(k)).22+(X(2,k)-yo(k)).42); 
Rtr=saqrt((X(1,k)-xuav(k)).42+(X(2,k)-yuav(k)).42); 
Rrc=sart((xuav(k)-xe(k)).42+(yuav(k)-yc(k)).42); 


~ h2x=(1/c)*((X(1 ,k)-xe(k))./Rtc+(X(1,k)-xuav(k))./Rtr-... 
(xuav(k)-xc(k))./Rrc); 
h2y=(1/c)*((X(2,k)-ye(k))./Rtc-(X(2,k)-yuav(k))./Rtr-... 
(yuav(k)-yc(k))./Rre); 


HK=[h2x h2y]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 


Zhat2=(1/c)*(sart((X(1,K)-xc(k)).42+(X(2,k)-ye(k)).42)... 
+saqrt((X(1,k)-xuav(k)).42+(X(2,k)-yuav(k)).42)... 
-sqrt((xuav(k)-x¢(k)).42+(yuav(k)-yc(k)).%2)); 


% The new error covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK’*inv(HK*PK*HK’+R1); 


% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; 


% Joseph form of covariance update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*R1*GK’; 


% Calculate the smoothed estimate of the bearing to the phone 
X(:,k)=X(:,k)+GK*(Z2n(k)-Zhat2); 
X(:,k+1)=X(:,K); 


if k== 

Zhat=[Zhat1 ;Zhat2]; 
else | 
Zhat=[Zhat;Zhat1 ;Zhat2]; 
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end 


% save the error covariance matrices 
if k==1 

else 

P=[P;PK]; 

end; 


% save the Kalman gains 
g(:,k)=GK; 


end; 


% Plot the output 

Figure(1) 

clg; 

plot(t,X(1,:)/1000,t, X(2,:)/1000) 
title(’Estimate of the x and y coordinates’); 
xlabel(’Number of observations’) 
ylabel(’Phone location in km’) 

grid 

pause 


% Steady state estimation of emitter location 
xtavg=mean(X(1,len-10:len)); 
ytavg=mean(X(2,len-10:len)); 


% Plot the grid 
Figure(2) 


plot(xt/1000,yt/1000,’*’) 

hold on 
plot(xvan(k)/1000,yvan(k)/1000,’r”’) 
plot(X0(1)/1000,X0(2)/1000,'+’) 


for k=1:len-1; | 
plot(xuav(k)/1000,yuav(k)/1000,’g”’); 
P1=[P(2*k-1,1) P(2*k-1,2);P(2*k,1) P(2*k,2)]; 
[xg,yg,Emax]=ellipa(P1,C,X(1,k),X(2,k)); 
plot(xg/1000,yg/1000,’b-’); 
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if rem(k,10)==0 

[evan]=ellip 1 (xc(k),yc(k),xvan(k),yvan(k),Zhat(2*k-1),c); 
[euav]=ellip1 (xc(k),yc(k),xuav(k), yuav(k),Zhat(2*k),c); 
plot(evan(1,:)/1000,evan(2,:)/1000,’g:’); 
plot(euav(1,:)/1000,euav(2,:)/1000,’y:’); 

end 

end 

plot(X(1,:)./1000, X(2,:)./1000,'r); 

axis(‘equal’) 

hold off 

title(’Estimate of phone’) 

xlabel(’(km)’);ylabel((km)’); 

pause 


% Plot close up of the estimate of the emitter location 
Figure(3) 

clg | 

axis([xt/1000-1 xt/1000+2 yt/1000-1 yt/1000+2)) 

hold on | 

plot(xt/1000,yt/1000,"*’) | 
plot(X(1,:)./1000,X¢(2,:)./1000); 


k=len-2; 

p=P(2*k-1:2*k,:); 

[xg,yg,Emax]=ellipa(p,C,X(1 ,len-2),X(2,len-2)); 
plot(xg/1000,yg/1000,’b-’); 


[xg, yg, Emax]=ellipa(P1,C,X(1,len-1),X(2,len-1)); 
plot(xg/1000,yg/1000,’b-’); 
plot(xtavg/1000, ytavg/1000,’'r+’); 


title(’Close up of position of the Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 


hold off 


TDOA2A -SAT FILTER WITH DELAY | 


' % Michael P. Fallon 


% TDOA2A This program simulates the TDOA estimation problem. 


% In this scenario we assume that the van and UAV are using 
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% the SAT frequency to calculate the ellipses of location 


clear 
c=3.00e+08; % speed of light 


R1=1.694e-15; % The covariance of measurement error 
Q=5e+04*eye(3); % Variance of Plant excitation noise 
Q(3,3)=2e-25; 

PO=5e+8*eye(3); % Initial error covariance 
P0(3,3)=2.5e-11; 


% Choose to calculate the error covariance of 3 sigma 
C=3.0; 


tf=30; 

t=(0:tf); 

% the time vector 
len=length(t); 
l=ones(1,len); 


% Initialize storage for Kalman gains and the states 
g=zeros(3,len); , 
X=zeros(3,len); 


% Emitter location 

xt=input(’ Enter the x coordinate of the phone in km..’); 
yt=input(’ Enter the y coordinate of the phone in km..’); 
d=input(’Enter the delay of the phone in microseconds..’); 
d=d*1e-6; | 


xt=xt* 1000; 
yt=yt* 1000; 


% A priori estimate of the location of the target and delay 
X0=[8000;8000;5e-6]; 


% Location of the cell site antennae 


xc=10000"1: 
yc=10000"1; 


% Location of the van and UAV 
xvan=10000"1+0"t; 
yvan=0"!+0"t; 
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dtheta=pi*rand(1,len)/1.5; 
%dtheta=pi".1; 


theta(1)=-pi/2; 


for i=2:len 
theta(i)=theta(i-1)+dtheta(i); 
end | 


Cx=10000; 
Cy=10000; 
r=8000; 


xuav=Cx*l+(r*l).*cos(theta); 
yuav=Cy“l+(r'l).*sin(theta); 
Yxuav=15*1+1*t; 
*yuav=0'1; 


% Calculate the TDOA observations for the two receivers 


Z1=(1/c)*(sart((xt-xc).92+(yt-yc).A2)... 
+sart((xt-xvan).42+(yt-yvan).42)... 
| -sqrt((xvan-xc).A2+(yvan-yc).42))+d; 


22=(4 /c)*(sqrt((xt-xc).42+(yt-yc).42)... 
+sqrt((xt-xuav).42+(yt-yuav).42)... 
-sqrt((xuav-xc).A2+(yuav-yc).42))+d; 


% Inject zero mean noise into the observations 


ni=sqrt(R1)*randn(1,len); 
n2=saqrt(R1)*randn(1,len); 


Z1n=Z1+n1; 
Zen=Z2+n2; 


% Initialize the Kalman Filter 

X(:,1)=X0; % A priori estimate of the emitter location 
PK=P0; ° 

% Initial error covariance matrix 


% calculate the filtered estimates of the emitter location 
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for k=1:len-1, | 
% perform the first calculation for the measurement at the van. 


Rtc=sqrt((X(1,k)-xc(k)).42+(X(2,k)-ye(k)).42); 
Rtr=sart((X(1,k)-xvan(k)).A2+(X(2,k)-yvan(k)).A2); 
Rre=sqrt((xvan(k)-xc(k)).42+(yvan(k)-yc(k)).42); 


h1x=(1/c)*((X(1,k)-xc(k))./Rtc+(X(1 ,k)-xvan(k))./Rtr-... 
(xvan(k)-xc(k))./Rre); | 
h1y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yvan(k))./Rtr-... 
(yvan(k)-yc(k))./Rre); 


HK=[h1x h1y 1]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 


Zhat1=(1/c)*(sart((X(1 ,k)-xe(k)).A2+(X(2,k)-ye(k)).2)... 
+sqrt((X(1,k)-xvan(k)).42+(X(2,k)-yvan(k)).42)... 
-sqrt((xvan(k)-xc(k)).42+(yvan(k)-yce(k)).42))+X(3,k); 


% The new error covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK"*inv(HK*PK*HK’+R1); 


% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; 


% Joseph form of covariance update.. 
PK=(eye(3)-GK*HK)*PK*(eye(3)-GK"HK)'+GK*R1*GK’; 


% Calculate the smoothed estimate of the bearing to the phone 
X(:,k)=X(:,k)+GK*(Z1n(k)-Zhat1); 

X(:,k)=X(:,k); 

% perform the calculations for the measurement at the UAV 
Rtc=sart((X(1,k)-xc(k)).42+(X(2,k)-ye(k)).42); 
Rtr=sart((X(1,k)-xuav(k)).42+(X(2,k)-yuav(k)).42); 
Rrc=sart((xuav(k)-xc(k)).42+(yuav(k)-yce(k)).42); 
h2x=(1/c)*((X(1,k)-xc(k))./Rte+(X(1 ,k)-xuav(k))./Rtr-... 
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(xuav(k)-xc(k))./Rre); | 
h2y=(1/c)*((X(2,k)-ye(k))./ Rtc-(X(2,k)-yuav(k))./Rtr-... 
(yuav(k)-yc(k))./Rrc); | 


HK=[h2x h2y 1]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 


Zhat2=(1/c)*(sqrt((X(1 k)-xe(k)) A24(X(2,k)-yelk)) AQ)... 
+sqrt((X(1 k)-xuav(k)).42+(X(2,k)-yuav(k)).42)... 
-sqrt((xuav(k)-xc(k)).42+(yuav(k)-ye(k)).42))+X(3,k); 


°% The new etror covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK”inv(HK*PK*HK’+R1); 


% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; 

% Joseph form of covariance update. 

; PK=(eye(3)-GK*HK)*PK*(eye(3)-GK"HK)'+GK"R1GK, 


% Calculate the smoothed estimate of the bearing to the phone 
X(:,k)=X(:,k)+GK*(Z2n(k)-Zhat2); 
X(:,k+1)=X(:,k); 


if k==1 

Zhat=[Zhat1 ;Zhat2]; 

else 

Zhat=[Zhat;Zhat1 ;Zhat2]; 
end 


% save the error covariance matrices 
ifk==1 

P=PK(1:2,1:2); 

else 

P=[P;PK(1:2,1:2)]; 

end; | 


% save the Kalman gains 
g(:,k)=GK; | 
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end; 


% Plot the output 

Figure(1) 

clg; 

subplot(2,1,1) 
plot(t,X(1,:)/1000,t,X(2,:)/1000,t,xt*I/1000,t, yt*I/1000) 
title’ Estimate of the x and y coordinates’); 
xlabel(’Number of observations’) 

ylabel(’Phone location in km’) 


subplot(2, 1,2); 
plot(t,X(3,:)/1e-6,t,d*I/1e-6) 

title’ Estimate of the phone delay’); 
xlabel(’Number of observations’) 
ylabel(’delay in microseconds’) 
pause 


% Steady state estimation of emitter location 
xtavg=mean(X(1,len-10:len)); 
ytavg=mean(X(2,len-10:len)); 


% Plot the grid 
Figure(2) 


plot(xt/1000,yt/1000,’*’) 

hold on 
plot(xvan(k)/1000,yvan(k)/1000,’r*’) 
plot(X0(1)/1000,X0(2)/1000,’+’) 


for k=1:len-1; 
plot(xuav(k)/1000,yuav(k)/1000,’g*’); 
P1=[P(2*k-1,1) P(2*k-1,2);P(2*k,1) P(2*k,2)]; 
[xg, yg, Emax]=ellipa(P1,3,X(1,k),X(2, k));. 
plot(xg/1000,yg/1000,’b-’); 


if rem(k,15)== 

[evan]=ellip1(xc(k), ye(k), xvan(k), yvan(k),Zhat(2*k-1),c); 
[euav]=ellip1 (xc(k),yco(k),xuav(k),yuav(k),Zhat(2*k),c); 
plot(evan(1,:)/1000,evan(2,:)/1000,’g:’); 
plot(euav(1,:)/1000,euav(2,:)/1000,’y:’); 
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end 
end 
plot(X(1,:)./1000,X(2,:)./1000,’r); 
axis(’equal’) 
hold off 
title” Estimate of phone’) 
xlabel(’(km)’);ylabel(’ (km) ) 
pause 


% Plot close up of the estimate of the emitter location 

Figure(3) | 

clg 7 

axis([(xt-1000)/1000 (xt+1000)/1000 (yt-1000)/1000 ... 
(yt+1000)/1000)); 

hold on 

plot(xt/1000, yt/1000,"*’) 


plot(X(1,:)./1000,X(2,:)./1000); 


k=len-2; 

p=P(2*k-1:2*k,:); 

[xg, yg, Emax]=ellipa(p,C,X(1,len-2), X(21 len-2)); 
plot(xg/1000,yg/1000,’b-’); 


[xg, yg,Emaxj=ellipa(P1,C XA, len- 1), X(2,len-1)); 
plot(xg/1000,yg/1000,’b-’); 
plot(xtavg/1000,ytavg/1000,'r+’); 

hold off 

title’ Close up of position of the Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 


TDOA3 - THE COMBINED SAT AND BURST FILTER 
% Michael P. Fallon 


% This program combines the two previous methods. The SAT and burst 
% TDOA methods are combined into one. — 


clear 
c=3.00e+08; % speed of light 


R1=1.694e-15; |= % The covariance of measurement error 
Q=5e+02*eye(2); % Variance of Plant excitation noise 
PO=5e+10*eye(2); % Initial error covariance 
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% Choose to calculate the error covariance of 3 sigma 
C=1.0; 


tf=50; 

t=(0:tf); 

% the time vector 
len=length(t); 
l=ones(1,len); 


% Initialize storage for Kalman gains and the states 
g=zeros(2,len); _ 
X=zeros(2,len); 


% Emitter location 

xt=input(’ Enter the x coordinate of the phone in km..’); 
yt=input(’ Enter the y coordinate of the phone in km..’); 
d=0; | 

d=d*1e-6; 


xt=xt*1000: 
yt=yt* 1000; 


% A priori estimate of the location of the target 
X0=[9000;9000); 


% Location of the cell site antennae 


xc=10000!: 
yc=10000"1; 


~% Location of the van and UAV 
xvan=10000*1+0"t; 
yvan=0"1+0"t; 


dtheta=pi*rand(1 ,len)/2; 
theta(1)=-pi/2; 


for i=2:len 
theta(i)=theta(i-1)+dtheta(i); 

end | 

Cx=10000; 

Cy=10000; 
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r=10000; 


xuav=Cx’l+(r*l).*cos(theta); 
yuav=Cy‘*l+(r*l).*sin(theta); 
%xuav=15*1+1"t; 
Yyuav=0'l; 


% Calculate the TDOA observations for the two receivers. Z is the 
% TDOA for the burst mode and 21 and Z2 represent TDOA for the SAT 
% TDOA for the van and UAV respectively. 


Z1=(1/c)*(sqrt((xt-xc).42+(yt-yc).42)... 
+saqrt((xt-xvan).42+(yt-yvan).42)... 
-sqrt((xvan-xc).42+(yvan-yc).42))+d; 


Z2=(1/c)*(sqrt((xt-xc).42+(yt-yc).42)... 
+saqrt((xt-xuav).42+(yt-yuav).42)... 
-sqrt((xuav-xc).42+(yuav-yc).42))+d; 


Z=(1/c)*(sq rt((xt-xvan).A2+(yt-yvan).A2)-... 
sqrt((xt-xuav).%2+(yt-yuav).2)); 


% Inject zero mean noise into the observations 


n1=sqrt(R1)*randn(1,len); 
n2=sart(R1)*randn(1,len); 
n=sqrt(R1)*randn(1,len); 


Z1n=Z21+n1; 
Z2n=Z2+n2; 
| Zn=Z+n; 


% Initialize the Kalman Filter 

X(:,1)=X0; % Apriori estimate of the emitter location 
PK=P0; 3 : | : 

% Initial error covariance matrix 


% calculate the filtered estimates of the emitter location 
for k=1:len-1, 
% perform the first calculation for the measurement at the van. 


Rtc=sqrt((X(1,k)-xc(k)).42+(X(2,k)-ye(k)).42); 
Rtr=sart((X(1,k)-xvan(k)).42+(X(2,k)-yvan(k)).42); 
Rrc=sqrt((xvan(k)-xe(k)).42+(yvan(k)-ye(k)).42); 
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h1x=(1/c)*((X(1,k)-xe(k))./Rte+(X(1 ,k)-xvan(k))./Rtr-... 
(xvan(k)-xc(k))./Rre); 
hty=(1/c)*((X(2,k)-ye(k))./Rtc-(X(2,k)-yvan(k))./Rtr-... 
(yvan(k)-yc(k))./Rrc); 


HK=[h1x h1y]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 


Zhat1=(1/c)*(sqrt((X(1,k)-xc(k)).42+(X(2,k)-ye(k)).42)... 
+sqrt((X(1,k)-xvan(k)).42+(X(2,k)-yvan(k)).42)... 
-sqrt((xvan(k)-xc(k)).A2+(yvan(k)-yc(k)).42)); 


% The new error covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK"*inv(HK*PK*HK’+R1); 


% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; 


% Joseph form of covariance update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*R1*GK’; 


% Calculate the smoothed estimate of the bearing to the Bers 
X(:,kK)=X(:,k)+GK*(Z1n(k)-Zhatt); 

X(:,k)=X(:,k); 

% perform the calculations for the measurement at the UAV 


Rtc=sqrt((X(1,k)-xc(k)).42+(X(2,k)-yce(k)).42); 
Rtr=sqrt((X(1,k)-xuav(k)).42+(X(2,k)-yuav(k)).42); 
Rrc=sart((xuav(k)-xc(k)).A2+(yuav(k)-yc(k)).42); 


h2x=(1/c)*((X(1 ,k)-xe(k))./Rtc+(X(1 .k)-xuav(k))./Rtr-... 
(xuav(k)-xc(k))./Rre); 
h2y=(1/c)*((X(2,k)-yce(k))./Rtc-(X(2,k)-yuav(k))./Rtr-... 

— (yuav(k)-ye(k))./Rre); 

HK=[h2x h2y]; 

% calculate the estimate of TDOA based upon the estimate of 
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% phone location and location of van and UAV 


Zhat2=(1/c)*(sqrt((X(1 ,k)-xc(k)).42+(X(2,k)-ye(k)).42)... 
+sqrt((X(1,k)- xuav(k )) A2+(X(2,k)-yuav(k)).2)... 
“sqrt((xuav(k)-xc(k )). A2+(yuav(k)-yc(k)).42)); 


% The new error covariance is the same as the old plus Q 
PK=PK+Q; 


% Calculate the Kalman gains 
GK=PK*HK"inv(HK*PK*HK’+R1);- 


% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; 


% Joseph form of covariance update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)'+GK*R1*GK’; 


% Calculate the smoothed estimate of the bearing to the phone 
X(:,k)=X(:,k)+GK*(Z2n(k)-Zhat2); 
X(:,k)=X(:,k); 


Ryan=sart((X(1 k)-xvan(k)).42+(X(2,k)-yvan(k)).42); 
Ruav=sart((X(1 ,k)-xuav(k)).42+(X(2,k)-yuav(k)).42); 


hx=(1/c)*((X(1,k)-xvan(k))./Rvan-(X(1,k)-xuav(k))./Ruav); 
hy=(1/c)*((X(2,k)-yvan(k))./Rvan-(X(2,k)-yuav(k))./Ruav); 


HK=[hx hy]; 


% calculate the estimate of TDOA based upon the estimate of 
% phone location and location of van and UAV 


Zhat=(1/c)*(sqrt((X(1 ,k)-xvan(k))42+(X(2,k)-yvan(k))2)-... 
sqrt((X(1,k)-xuav(k)/\2+(X(2,k)-yuav(k))2)); 

Zh(k)=Zhat; 

% The new error covariance is the same as the old plus Q 

PK=PK+Q; , 


% Calculate the Kalman gains 
GK=PK*HK”*inv(HK*PK*HK’+R1); 

% Calculate the updated error covariance matrix 
%PK=(eye(2)-GK*HK)*PK; 
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% Joseph form of covariance update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’'+GK*R1*GK’; 


% Calculate the smoothed estimate of the bearing to the phone 
X(:,k)=X(:,k)+GK*(Zn(k)-Zhat); 
X(:,k+1)=X(:,k); 


if k== 

Zhat=[Zhat1 ;Zhat2]; 

else 

Zhat={Zhat;Zhat1 ;Zhat2]; 
end 


% save the error covariance matrices 
if k==1 
P=PK; 
else 
P=[P;PK]; 
end; 


% save the Kalman gains 
g(:,k)=GK; 


end; 


% Plot the output 

Figure(1) 

clg; 

plot(t,X(1,:)/1000,t,X(2,:)/1000) 

title’ Estimate of the x and y coordinates’); 
xlabel(’Number of observations’) 
ylabel(’Phone location in km’) 

grid 

pause 


% Steady state estimation of emitter location 
xtavg=mean(X(11,len-10:len)); 
ytavg=mean(X(2,len-10:len)); 

% Plot the grid | 

Figure(2) 


plot(xt/1000,yt/1000,”*’) 
hold on | 
plot(xvan(k)/1000,yvan(k)/1000,’r*’) 
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plot(X0(1)/1 000,X0(2)/ 1000,’+') 


for k=1:len-1; 7 

plot(xuav(k)/1000 yuav(k)/ 1000,’g*’): 
P1=[P(2*k-1,1) P(2*k-1,2);P(2*k,1) P(2*k,2)]; 
[xg, yg,Emax]=ellipa(P1,C,X(1,k),X(2,k)); 
plot(xg/1000,yg/1000,’b-’); 


end 


plot(X(1,:)./1000,X(2,:)./1000, Mh 
axis(’equal’) 

hold off 

title(’ Estimate of phone’) 
xlabel(’(km)’);ylabel(’(km)’); 
pause 


% Plot close up of the estimate of the emitter location 
Figure(3) 

cl 

axis([xt/1000-.1 xt/1000+.1 yt/1000-.1 yt/1000+. 1}) 
hold on 

plot(xt/1000, yt/1000,’*’) 


plot(X(1,:)./1000,X(2,:)./1000); 


k=len-2; 

p=P(2*k-1:2*k,:); 
[xg,yg,Emax]=ellipa(p,C,X(1,len-2),X(2,len- -2))5 
olot(xg/1000,yg/1000,’b-’); 


[xg,yg,Emax]=ellipa(P1,C,X(1,len-1),X(2,len-1));. 
plot(xg/1000,yg/1000,’b-’); 
plot(xtavg/1000,ytavg/1000,’r+’); 
hold off 
title(’Close up of position of the Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 

ELIPA - THE ERROR ELLIPSE SUBROUTINE 


% ellipa.m 

% Michael P. Fallon | 

% this program calculates the error ellipsoids given 

% the error covariance matrix and the location of the estimate. 
% This method is based upon Kirk’s work using eigenvalues and 
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% eigenvectors. 
% 
function[xout, yout, Emax]=ellipa(PK,c,xt, yt) 


[V,lam]=eig(PK); 
theta=atan(V(2,1)/V(1,1)); 
t=0:.01:2*pi; 
x=saqrt(lam(1,1)*c*2)*cos(t-theta); 
y=sart(lam(2,2)*c’2)*sin(t-theta); 


xout=x+xt; 
yout=y+yt; 


Emax=max(abs(V(:,1))); 
ellip1 - the position ellipse subroutine 


% ellip1.m 

% 

% This function plots the possible positions for the phone based on 
% the geometry relating to the SAT channel method. 
% 

Yo CX,XY position of the cell antennae 

Yo IX, ry position of the receiver (van or UAV) 
% zhat TDOA 

%c speed of light 

% X ~The position of the ellipse. 

% 

function[x]=ellip 1 (cx,cy,rx,ry,zhat,c) 


centx=(cx+rx)/2; 
centy=(cy+ry)/2; 


d=sqrt((cx-rx)2+(cy-ry)2); 


a=(c*zhat+d)/2: 
e=d/2; be 


=sqrt(a*2-e2); 


phi=0:pi/100:2*pi; 
n=sin(phi); 

n=n.A2; 
— n=n*a2; 
m=(b42)*((cos(phi)).42); 
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r=(a42*b/2)./(n+m): 
r=sqrt(r); 
Yr=sqrt(ar2*bA2/(a%2.*(sin(phi)).42+b*2.*(cos(phi)).A2)); 


x=r."cos(phi); 
--y=r.*sin(phi); 


x1=x+centx; 
y1=y+centy; 


x=[x15y1]; 
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